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ABSTRACT. This article presents uniform random generators of plane partitions 
according to the size (the number of cubes in the 3D interpretation). Com- 
bining a bijection of Pak with the method of Boltzmann sampling, we obtain 
Q\ random samplers that are slightly superlinear: the complexity is 0(n(lnn) 3 ) 

C ~ ) , in approximate-size sampling and 0(ra 4 / 3 ) in exact-size sampling (under a 

C ~ ) ■ real-arithmetic computation model). To our knowledge, these are the first 

£N| ' polynomial-time samplers for plane partitions according to the size (there ex- 

ist polynomial-time samplers of another type, which draw plane partitions 
that fit inside a fixed bounding box). The same principles yield efficient sam- 
plers for (a X 6)-boxed plane partitions (plane partitions with two dimensions 
bounded), and for skew plane partitions. The random samplers allow us to 
perform simulations and observe limit shapes and frozen boundaries, which 
have been analysed recently by Cerf and Kenyon for plane partitions, and by 
Okounkov and Reshetikhin for skew plane partitions. 
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Plane partitions, originally introduced by A. Young [30j . constitute a natural 
generalisation of integer partitions in the plane, as they consist of a matrix of inte- 
gers that are non-increasing in both dimensions (whereas an integer partition is an 
array of non-increasing integers). In addition, they also have a nice interpretation 
■ in 3D-space as a heap of cubes (see Figure [2]) ■ Plane partitions have motivated a 

huge literature in numerous fields of mathematics [IJ [Til H21 El ESj and statistical 
physics [16] [27] , and have provided crucial insight for solving challenging problems 
in combinatorics [51], see [2| for a detailed historical account. The problem of enu- 
merating plane partitions was solved by MacMahon [T5] , who proved the beautiful 
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formula 



for the generating function. The simplicity of the formula asks for a combinatorial 
interpretation. A first direct bijective proof has been given by Krattenthaler [15] . 
The principle is inspired by the seminal bijection of Novelli-Pak-Stoyanovskii [5JJ 
giving an interpretation of the hook- length formula. In [T5J, Krattenthaler also 
discusses as application of his bijection a polynomial-time algorithm for the random 
generation of plane partitions in a given a x b x c box. Upon looking at the heap 
of cubes in the (1,1,1) direction, this task is equivalent to sampling tilings of a 
hexagon of side lengths (a, 6, c, a, 6, c) by rhombi; there also exist random samplers 
for such tilings, which rely either on "coupling from the past principles" [26j or 
on "determinant algorithms" |28| . In contrast, we are interested here in sampling 
plane partitions uniformly at random with respect to the size, defined as the sum 
of the matrix entries. For this purpose, we use another bijective interpretation of 
MacMahon's formula recently given by Pak [24 . 

Let us briefly mention the motivations for having a random sampler of plane par- 
titions according to the size. The size is a natural parameter, as it corresponds to 
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the volume of the plane partition (number of cubes) in the 3D interpretation. Re- 
cently, several authors have studied the statistical properties of plane partitions with 
respect to the size. In particular, under fixed-size distribution, Mutafchiev |18j has 
shown a limit law for the maximal entry, and Cerf and Kenyon [3] have determined 
the asymptotic shape (the asymptotic shape in the boxed framework — hexagon 
tilings — is due to Cohn, Larsen, and Propp [1]). Even more recently, Okounkov 
and Rcshctikhin, using a method based on Schur processes, have rediscovered the 
limit shape of Cerf and Kenyon [23] . They have studied in a subsequent article [53"] 
the local correlations and limit shapes for plane partitions under a mixed model: 
the plane partition is constrained to a 2-coordinate a x b box and is drawn under the 
Boltzmann model with respect to the size. (We will also describe random samplers 
for this mixed model.) In addition, physicists have developed new models relying 
on plane partitions, giving rise to a simplified version of the 3-dimensional models 
of lattice vesicles [14] . Plane partitions are also related to the 3-dimensional Ising 
model in the cubic lattice [3J. In general, physicists are interested in checking ex- 
perimentally or conjecturing some limit properties of these models, by generating 
very large random objects. 

For this purpose, this paper introduces efficient samplers for plane partitions. 
Our approach combines methods from bijective combinatorics and symbolic com- 
binatorics. Precisely, a minor reformulation of Pak' bijection maps a multiset of 
integer pairs (the class is denoted by A4) to a plane partition with the same size. 
(Since the class M. has generating function rir>i(^ — ' T ' )~ r ' gives a direct proof 
of ([T]).) Our aim here is to take advantage of this bijection for random sampling. 
Indeed Pak's bijection reduces the task of finding a sampler for plane partitions 
to the task of finding a sampler for M.. As the class M has an explicit simple 
combinatorial decomposition, it is amenable to random sampling methods from 
symbolic combinatorics. By now there is the recursive method [20, 9\ based on the 
counting sequences and Boltzmann samplers based on the generating functions, 
as introduced in [5] and further developed in [BJ. We adopt here the framework 
of Boltzmann samplers, which tend to be more efficient as they avoid the costly 
precomputations of coefficients required by the recursive method. 

As opposed to the recursive method — which produces exact-size samplers — 
the probability distribution in Boltzmann sampling is spread over the whole class; 
precisely an object of size n has probability proportional to x n , where x is a fixed 
real parameter. In particular, as two objects having the same size have equal 
probability, the probability distribution restricted to a given size n is uniform. As we 
are interested in generating very large plane partitions, the Boltzmann framework 
is suitable, due to the gain obtained by relaxing the exact-size constraint. The 
articles [5j and [BJ provide a collection of rules for building a sampler for a class 
admitting a decomposition involving classical constructions. Using these rules, the 
decomposition of M. is readily translated into a Boltzmann sampler. This yields, 
via Pak's bijection, a Boltzmann sampler for plane partitions. In addition, as 
the size distribution of plane partitions — under the Boltzmann model — has good 
concentration properties, it is possible to "tune" the parameter x so as to draw 
objects of size around (or exactly at) a given target value n. With the parameter 
x suitably tuned and a rejection loop targeted at the size, we obtain a quasi- 
linear time approximate-size sampler for plane partitions: for any tolerance-ratio 
e £ (0, 1), our sampler draws a plane partition of size in [n(l — e), n(l + e)] with an 
expected running 0(n(ln n) 3 ). The same principles, with the rejection loop running 
until a given size n is attained, yields an exact-size sampler for plane partitions, 
with expected running time 0(n 4 / 3 ). To our knowledge, our algorithm is the first 
exact-size sampler for plane partitions with expected polynomial running time. This 



RANDOM SAMPLING OF PLANE PARTITIONS 



3 



allows us to generate objects of size up to 10 7 in a few minutes on a PC. The same 
principles (i.e., Pak's bijection + Boltzmann samplers) yield efficient Boltzmann 
samplers for (a x 6)-boxed plane partitions (plane partitions whose non-zero entries 
lie in an (a x b) rectangle) , which are those considered by Okounkov and Reshetikhin. 
We obtain for boxed plane partitions an approximate-size sampler with expected 
running time O aj b ]E (l) and an exact-size sampler of expected running time O a ^{n), 
where e is a tolerance-ratio on the size (for approximate-size sampling) and where 
n is the target-size (the dependency in a, b, e of the asymptotic constants in the big 
O's are stated precisely in Theorem [T0|) . 

Proving the correct complexity orders of the samplers is the major technical 
difficulty we have to deal with. At first we have to analyse the expected running 
time of the Boltzmann sampler for the multiset-class Ai, as well as the size dis- 
tribution on A4 under the Boltzmann model. All this is done using the Mellin 
transform. Second, we study the complexity of Pak's bijection, which depends on a 
natural length-parameter of a plane partition, which is the maximum hook-length 
(abscissa+ordinate+1) over all nonzero entries of the matrix. Let us finally mention 
that, for the sake of simplicity, all complexity results are stated and proved with 
the O notation (upper bound). With little more care one could prove that all the 
stated complexity results hold in fact with a notation, i.e., an upper and a lower 
bound. 

Outline of the paper. After some definitions in Section [1] about combinatorial 
classes and plane partitions, we present in Section [2] a slight reformulation (more 
algorithmic) of Pak's bijection. The bijection induces a combinatorial isomorphism 
between the set of plane partitions and the class M := MSet(Z x Seq(Z) 2 ). 
Section [3] recalls basic principles of Boltzmann sampling, in particular the sampling 
rules associated to the constructions appearing in the specification of A4. The 
Boltzmann sampler for M, as well as M a .b '■= MSet(Z x Seq <q (Z) x SEQ <h (Z)), 
is derived in Section [21 giving rise to Boltzmann samplers for plane partitions and 
(a x 6)-boxed plane partitions. We explain then briefly how the principles extend 
to obtain Boltzmann samplers for so-called skew plane partitions. In Section l4~4l 
using suitable choices of the parameter x in the Boltzmann samplers, we obtain 
efficient samplers for plane partitions targeted exactly or approximately at a given 
size n (precise statements are given in Theorems I§1 and 1 10|) . The expected running 
times of the targeted samplers are then analysed in Section [H 

1. Definitions 

A combinatorial class is a pair (A, |.|) where A is a set and |.| is a function from 
A to N, called the size function, such that the number of elements of any given 
size is finite. Using the size function, we can graduate A as A — U n A n , where 
A n is the set of objects of A that have size n. In the sequel, we denote by A n 
the cardinality of A n - To each combinatorial class A, we associate the generating 
function A{z) = ^2A n z n . 

Two combinatorial classes (A, and (B, |.|g) are said to be combinatorially 
isomorphic, A ^ B, if and only if there exists a one-to-one mapping from A to B 
that preserves the size. Let us notice that two classes A and B are isomorphic if 
and only if their generating functions are equal. 

Here are some classical constructions on combinatorial classes that will be used 
in this paper. Notations and rules are summarized in Figure [T] (a more general 
presentation can be found in [8]): 

- £ and Z are atoms of size and I. 

- Disjoint union A + B: the union of two copies of A and B made disjoint. 

- Cartesian product A x B: the set of pairs (a, (3) where a £ A and (3 e B. 
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Given a class A not containing empty atoms, 

- Sequence: Seq(^4) is the class of finite sequences of objects of A. 

- Multiset: MSet(^) is the class of finite sets of objects of A, with repetitions 
allowed. 

In all these constructions, the size of an object in the composed class is naturally 
defined as the sum of the sizes of the components (e.g., the size of a sequence 
71, . . . ,7fc is |7i| + • • • + |7fc|)- Observe that, in a multiset /x € MSet(^4), each 
element a £ A has a multiplicity c a > 0. Hence, if A is a finite set, 

(2) MSet(^) ~ Y[ SEQ({a}). 

aeA 



Class 


Generating function 


Definition 


C = £ 
C = Z 


C(z) = 1 
C{z) = z 


neutral object of size 
atom of size 1 


C = A + B 
C = AxB 
C = Seq(^) 
C = MSet(„4) 


C(z) = A(z) + B{z) 
C{z) = A(z) x B(z) 
C(z) = (l-A(z))- 1 
C(z) = exp(EMz k )/k) 


disjoint union 

cartesian product 

£ + A + AxA + AxAxA + ... 

a multiset of elements of A 



Figure 1 . Some constructions on combinatorial classes. 

A plane partition (Figure H|) of n is a two-dimensional array of non-negative 
integers (aij)^2 that are non-increasing both from left to right and bottom to top 
and that add up to n. In other words, 

(3) V(i,j) € N 2 a it j > Oij+i, a it j > a i+ ij and a i:j = n. 

We denote by V the combinatorial class of plane partitions, endowed with the size 
function |(aij) N2 | = ^2ij a i,j- Plane partitions have a natural representation in 
3D-space as a heap of cubes with non-increasing height in the direction of the x- 
axis and y-axis, see Figure [H Observe that the size of the plane partition exactly 
corresponds to the number of cubes in the 3D-representation. 

The bounding rectangle of a plane partition (a* j )n 2 is the smallest double range 
R = [0.1 — 1] x [Q..w — 1] such that dij = for all index pairs outside of R. 
An (a x b) -boxed plane partition is a plane partition whose bounding rectangle is at 
most a x b. Equivalently, djj is null for any such that i > a or j > b. We denote 
by V a ,b the class of (a x 6)-boxed plane partitions. Define the two combinatorial 



Figure 2. Plane partition of size 22 and its 3D representation. 

classes M and M a ,b as follows, where SEQ <d (^4) denotes the class of sequences of 
at most d — 1 elements of A. 

(4) M := MSet(Z x Seq(Z) 2 ) 

(5) M a .b := MSet(Z x Seq <q (Z) x SEQ <b (Z)). 
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Classically, Seq(Z) is identified with the class of nonnegative integers, so that 
we can specify M. with the following simplified notation, 

(6) M ~ MSet(Z x N 2 ). 

In the next section we explain how Pak's bijection yields an explicit combinatorial 
isomorphism between V and M. For this purpose, we introduce some more termi- 
nology. The diagram D of an element M 6 M. is a two-dimensional array (mjj)^ 
(with (0, 0) at the bottom left) where m^j is the multiplicity of (Z, i, j) in M (see the 
first two pictures of Figured!]). The size of D is defined as \D\ = J2i j TO ij(*+.? + l)i 
so that it corresponds to the size of the multiset in A4. The bounding rectangle 
of M is defined similarly as for plane partitions: it is the smallest double range 
R = [0.4 — 1] x [0..w — 1] such that all entries of D outside of R are zero. The 
integers £ and w are respectively called the length and the width of D. Note that, 
for fixed integers a and b, the diagrams of elements in M. a ,b are constrained to have 
their bounding rectangle C [0..a — 1] x [0..6 — 1]. Therefore M. a ,b is called the class 
of (a x b)-boxed multisets. 



2. Pak's bijection 

In |24j . Pak presents a bijection between plane partitions bounded in a shape fj, (/x 
being a Ferrers diagram) and fillings of the entries of /j, with nonnegative integers. 
We reformulate this bijection as an algorithm, Algorithm [T] below, that realises 
explicitly the combinatorial isomorphism M. ~ V, see Figure |3] for an example. 



Algorithm 1: From the diagram of a multiset to a plane partition 
Input : The diagram D of a multiset in M. 
Output: a plane partition. 
Let I be the length and w be the width of D. 
for i <— I — 1 downto do 

for j ' ^- w — 1 downto do 

D[i,j] «- D[i,j] + ma X (D[ l + l,j]),D[i,j + 1]); 
for c <— 1 to min(™ — 1 — i, £ — 1 — j) do 
x «- i + c; y «- j + c; 
|_ max(D[a;+l, y], D[a;, y+l])+ min(D[s-l, y], D[a;, y-l])-£)[as, y]; 

return D; 



Proposition 1 (Pak [24j ). Algorithm]^ yields an explicit size-preserving bijection 
between the class fA a ,b cmd the class of (a x b)-boxed plane partitions. In other 
words, the algorithm realises the combinatorial isomorphism 

(7) V a , h ~ MSET(Z X SEQ <a (Z) X SEQ <b (Z)). 

Proof. See [24]. □ 

Proposition 2. Algorithm]]] realises the combinatorial isomorphism 

(8) V ~ MSet(Z x Seq(Z) 2 ). 

Proof. Take the limit a — > oo and & — > cxd in Proposition [T] □ 

3. Boltzmann sampling 

This section recalls basic principles of approximate-size sampling under Boltz- 
mann model (0 [(>])• 
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Figure 3. Pak's bijection on an example. 



Definition 3 (Boltzmann model). Let C be a combinatorial class and C(x) :— 
J-^ec^' 7 ' its generating function. Given a coherent positive real value of x, i.e., 
chosen within the disk of convergence of C{x) ), the Boltzmann model of parameter 
x assigns to any element 7 S C the following probability, 

rlfl 

A Boltzmann sampler TC(x) for C is an algorithm that produces objects of C 
at random under the Boltzmann model. As elements of the same size have the 
same weight, the probability induced by a Boltzmann sampler on any given size n 
is uniform. The size of the output is a random variable N x satisfying 

C n x n 



C(x) 



Figure 0] shows this probability distribution for plane partitions. When a target 
size n has to be achieved, the idea is to tune the parameter x so that M(N X ) = n 
(see Section [5]). 



P(7V r =«) 




FIGURE 4. Probability distribution of sizes for plane partitions under 
Boltzmann model, with different values of the parameter x. 



Figure [5] briefly summarizes how to obtain samplers for the constructions used 
in the specification of M. (see details in [5]); the rules can be combined to build 
a generator for any class specified with these constructions, in particular the class 
A4. The sampling rules make use of simple auxiliary generators: Geom(p) generates 
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C = A x B 


FC(x) := {TA(x),FB(x)} 


C = Seq(^) 


TC(x) := [Geom(A(x)) =^ TA(x)] 

where G FY means " return G independent calls to FY" . 


C = MSet(^) 


Define the probability distribution relative to A and x: 

Pr(K < k) = Y\ exp ( - -^ar*)). 

]>k 3 

Let MAX_lNDEx(yl; x) be a generator according to this distribution. 

TC(x) : 7<— 0; ko <— Max_Index(A; x); 
if fco then 

for j from 1 to fco — 1 do 

p^Pois^f ); 

for i from 1 to p do 

7^7, copy(ry4(a; J ) j times) 

p^Pois^^l); 
for i from 1 top do 

7 ^— 7, copy(r J 4(a;' C0 ) fco times) 
return 7. 



Figure 5. Sampling rules associated to Boltzmann samplers for 
some combinatorial constructions. 



integers under the geometric law P(fc) = p k (l — p), Pois(A) generates integers under 
the Poisson law P(fc) = e _A ^j, and Pois>i(A) generates integers under the positive 

Poisson law P(fc) = (e A — ^)~ 1 jt for k > 0. Such generators are classically realised 
by simple iterative loops, the complexity of generating an integer k being 0(k), 
see [5] for a discussion. 

Proposition 4 (Flajolet et al. [5]). Given two combinatorial classes A and B 
endowed with Boltzmann samplers T A{x) andTB{x), the sampler TC(x), as defined 
in the first entry of Figure is a Boltzmann sampler for A x B. Given a class 
A not containing the empty atom and endowed with a Boltzmann sampler TA(x), 
the samplers TC(x), as defined in the second and third entr$\ of Figure are 
respectively Boltzmann samplers for Seq(^4) and for MSet(A). 

From these sampling rules, a class C recursively specified from atomic sets in 
terms of these constructions can be endowed with a Boltzmann sampler TC(x). 
The complexity of generating an object 7 £ C is 0(|7|)- 

Complexity model. Let us say a few words on the specific real-arithmetic com- 
plexity model used for Boltzmann samplers. First, notice that the samplers given in 
Figure [5] draw integers according to distributions (Geom, Pois, MaxJndex) that 
require the exact values of the generating functions of the classes involved. Hence, 
such generating functions should be evaluated. The complexity model we adopt, 
as already defined in [5], relies on the oracle assumption. This assumption allows 
us to separate the combinatorial complexity of the sampler from the complexity of 
evaluating the generating functions (there are already some results [2 5) and work 
in progress dedicated to the latter issue). Given any combinatorial class C specified 
recursively using the constructions of Figure El and given a value x > within 
the disk of convergence of C(x), we assume that an oracle provides, at unit cost, 
the exact values at x of the generating functions for all classes intervening in the 
decomposition of C. In practice, we work with a fixed precision (e.g., 20 digits) 



^We hereby correct an omission in the definition of the sampler for MSET(yt) given in [6], 
namely the test fco 7^ 0. 
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and precomputc the values of generating functions used by the Boltzmann sampler. 
Finally let us come back to the complexity of drawing an integer K under a certain 
discrete distribution 

F(K = k)= p k , 

where the constant pk are known (if pk involves exact values of generating functions, 
the oracle provides the values of pk ) ■ As discussed in [5] , drawing an integer K under 
this distribution is done by a simple loop, 



Algorithm 2: Drawing an integer K under an arbitrary distribution 

U := uniform(0,l); S := 0; k : = 0; 
while U < S do 

L S:=S+ Pk ; k:=k + h 
return k; 



Thus, the cost of drawing K is of the order of the value k that is finally assigned 
to K . An exception is the case of the geometric law, which is simpler. Indeed, to 
draw K under Geom(a;) (with x £ [0,1]), it is enough to set K = [\n(U)/ ln(x)J , 
where U is uniform in (0, 1). 

Hence, the cost of drawing a geometric law is 0(1). 

4. Samplers for plane partitions 

4.1. Boltzmann sampler for plane partitions. The explicit bijection between 
M. and V allows us to design a simple Boltzmann sampler for plane partitions, 
made of two steps: (i) generate a multiset in A4 under the Boltzmann model, (ii) 
apply Algorithm [T] (Pak's bijection) to the diagram of the multiset generated. 



Algorithm 3: TM(x) [Boltzmann sampler for A4] 
M is the diagram of the multiset to be generated 
V{x,y) M[x,y] ^ 0; 

ko <— Max_Index(A; x), where A(x) = x/(l — x) 2 ; 
if ko 7^ then 

for k <— 1 to ko — 1 do 
p^Pois( fc(1 ^ fc) . ); 
for i *— 1 to p do 

x ^— Geom(x fc ); y «— Geom(x k ); 
|_ M[x,y] <- M[x,y] + k; 

p^ Pois >i( fc0 (i-°*o)2 ); 

for i *— 1 to p do 

x <- Geom{x k °); y <— Geom(x fco ); 
[ M[x,y]^ M[x,y] + k ; 

return M; 



Lemma 5. Given < x < 1, the generator TM(x) — as defined in Algorithm^ — 
is a Boltzmann sampler for M . 

Proof. The specification of M. , given in Equation ((6|) , is translated to a Boltzmann 
sampler using the rules of Figure [5j The translation is carried out directly on the 
diagram of the multiset (recall that the entry of the diagram corresponds to 
the multiplicity of (Z,i, j) in the multiset). □ 
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Algorithm 4: TP(x), with < x < 1 [Boltzmann sampler for plane partitions] 

Compute (i <— rM(ic); 

Apply Algorithm [1] (Pak's bijection) to n; 

return n 



Since Pak's bijection preserves the size, Algorithm 2] is a Boltzmann sampler for 
plane partitions. 

Figure E] shows computation time^l of TP(x) for sizes up to 10 7 : the first line 
gives the time of generation of the multiset (TM(x)) and the second line gives the 
computation time of Pak's bijection. The sampler has been implemented in Maple 
and the bijection in OCaml. As we can see, the complexity is dominated by Pak's 
bijection for objects of large size. This is confirmed by the analysis to be given in 
Section [SI the complexity of drawing a multiset of size around n is of order n 2/3 , 
while the expected running time of Pak's bijection applied to a random multiset of 
size n is of order n(lnn) 3 . Figure [7] shows two large plane partitions generated by 
TP(x) for x close to 1, x = 0.947 and x = 0.9866. 



approx. size 


10 3 


10 4 


10 5 


10 6 


10 7 


generation 


~ 0.1 sec. 


~ 0.5 sec. 


~ 2-3 sec. 


~ 10 sec. 


~ 60 sec. 


Pak's transform 


~ 0.1 sec. 


~ 0.3 sec. 


~ 2 sec. 


~ 20-30 sec. 


~ 8-9 min 



Figure 6. Time per generation for different sizes of plane partitions. 




(a) A random piano partition (b) A random plane partition of size 1,005,749 generated 
of size 15,256, generated by by rP(0.9866), seen from the direction (1, 1, 1). 

TP(0.947). 



Figure 7. 

4.2. Boltzmann sampler for (a x 6)-boxed plane partitions. According to 
the equivalence with the definition in terms of diagrams, an element of M. a ,b is a 



Computations have been performed on a Mac OS X Power PC G4 1,42GHz, with 1GB of 
RAM and 512 kB of cache. 
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multiset of pairs with < i < a and < j < b, each element having 

size (i + j + 1). The set of such pairs being finite, Equation (J2]) yields 

(9) M a . b = ;Q Seq(Z^' +1 ) 

0<i<a 
0<j<6 

Lemma 6. Given < x < 1, the generator TM a ^{x) — as defined in Algorithm^ — 
is a Boltzmann sampler for A4 a ,b- 

Proof. Translate the specification © to a Boltzmann sampler for M. a j> using the 
rules of Figure [5] □ 

Algorithm 5: rM a: j,(i) [Boltzmann sampler for M a ,b] 
M is the diagram of the multiset to be generated 
for i <— to a — 1 do 
for j *— to b — 1 do 

|_ M[i,j] <- Geom(2: i+j+1 ); 

return M; 

Again, since Pak's bijection preserves the size, the following generator is a Boltz- 
mann sampler for (a x 6)-boxed plane partitions. 

Algorithm 6: TP a ^(x), with < x < 1 [Boltzmann sampler for boxed plane 
partitions] 

Compute n <— TM ai b(x); 

Apply Algorithm [T] (Pak's bijection) to fi; 

return n 



4.3. Extension to skew plane partitions. We consider here a natural gener- 
alisation of (a x fe)-boxed plane partitions, called (a x 6)-boxed skew plane par- 
titions. A (a x 6)-boxed skew plane partition is given by an index-domain D C 
[0..a — 1] x [0..6 — 1] such that D is obtained from [0..a — 1] X [Q..b — 1] by removing 
rectangles of the form [0..a ! — 1] x [0..6'— 1], with a' < a and b' < b. Each truncation 
by a smaller rectangle makes an outer corner appear in the index domain (e.g., the 
partition of Figure [5] has 2 outer corners). Let us denote by Vd the class of all such 
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Figure 8. A skew plane partition of size 43. 

partitions for a given domain D. For this new class of partitions, we need to define 
the hook-length of a pair in the domain D. Let £(i) be the minimum abscissa 
such that £ D and d(j) the minimum ordinate such that (i, d(j)) G D. The 

hook-length of (i,j) in D is then h(i, j) = (i — £(i)) + (j — d(j)) + 1, which is exactly 
i + j + 1 when D is [0..a — 1] x [0..b — 1]. In [53], Pak's bijection is most gener- 
ally described for skew plane partitions, which leads to the following combinatorial 
isomorphism: 

Vn - } ; SEQ(Z h ^) 
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The Boltzmann sampler for (a x 6)-boxed plane partitions extends directly to a 
Boltzmann sampler for skew plane partitions as follows: to sample a diagram, 
draw the value at each point in D according to a geometric law of parameter 
x h(i,j) . t nen apply Algorithm Q] (Pak's bijection) to the multiset generated, with 
the difference that the domain scanned by is D. 

Okounkov and Reshetikhin J3] have studied the limit shape of a skew plane par- 
tition under the Boltzmann distribution, with the Boltzmann parameter x tending 
to 1. If the lengths of the rectangles are of order (1 — x) -1 , some interesting phe- 
nomena are to be observed regarding the typical shape of a random skew plane 
partition. Using a technique based on Schur processes, the authors of |23j provide 
a precise analysis of these phenomena. They prove that the (rescaled) limit shape 
of a skew plane partition has a frozen boundary that satisfies explicit equations, 
and they classify the non-smooth points of the boundary as turning points and 
cusps. Turning points always appear, even for a boxed domain (a x b); they cor- 
respond to points of tangency of the frozen boundary with the delimiting 3D-box 
{(x,y,z) G R+,s.i. (x,y) G D}. If the index domain has outer corners, some cusp 
points possibly appear at each of the outer corners. 

Our random sampler for skew plane partitions makes it possible to perform sim- 



ulations and observe these asymptotic phenomena. Figure 9(a) shows a (100 x 100)- 
boxed plane partition of size 999,400. A frozen boundary appears that meets the de- 
limiting 3D-box in a tangential way (these points of tangency are precisely the turn- 



ing points in the terminology of Okounkov and Reshetikhin). And Figure 9(b) shows 
a skew plane partition of size 1, 005, 532 on the index-domain (100 x 100)\(50 x 50), 
which has an outer corner at (50,50); accordingly a cusp point appears on the 
boundary of the limit shape at the point above (50, 50). 

Note that the typical shape of a large unconstrained random plane partition, as 
shown in Figure [7J has different features: there are 3 "legs" — one in each axis- 
direction — whose lengths tend to infinity even when the plane partition is rescaled 
to have unit volume (which essentially corresponds to rescaling by a factor (1 — x)^ 1 
in each dimension). 

4.4. Samplers targeted around a given size. Given a class C — U„C n , an 
exact-size sampler is a procedure that, for any given n > 1 (called the target-size), 
outputs an object of C n uniformly at random. An approximate- size sampler is a 
procedure that, for any given n > 1 and e G (0,1) (called the tolerance-ratio), 
outputs an object of C of size in [n(l — e),n(l + e)] and such that two objects of 
the same size have the same chance to be drawn (hence the distribution induced 
on each size k G [n(l — e), n(l + e)] is uniform). 

Such procedures are easily obtained if C is endowed with a Boltzmann sampler 
TC(x). Fix a suitable value of x and repeat calling TC(x) until the size is in the 
desired size-range fl; ft = {n} for exact-size sampling and 17 = [n(l— e), n(l+e)] for 
approximate-size sampling. The exact-size sampler and approximate-size sampler 
defined in this way are denoted SampleC(:e; n) and SampleC(x; n, e), respectively. 

In general one chooses x so that the expected size AC(x) of the output of TC(x) 
— which satisfies AC{x) = xC (x)/C(x) as proved in [5] — is equal to n, or at least 
is asymptotically equal to n, that is, one looks for an exact or an approximate 
solution of the so-called target- size equation: 

As we show in Section [5] the expected size AM{x) of the output of TM{x) satisfies 
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(a) A (100 X 100)-boxed plane parti- 
tion of size 999,400 drawn under 
Boltzmann distribution at X = 
0.9931. 




(b) A [0..99] x [0..99]\[0..49] x [0..49] 
skew plane partition of size 
1,005,532 drawn under Boltzmann 
distribution at x = 0.9942. 



Figure 9. 



where £(.) is the Riemann zeta function, so a suitable value of the parameter x to 
reach a target-size n is £ n = l — (^((fy/n) 1 / 3 , since 2£(3)/(l— £„) 3 = n; and we show 
in Section [5] that the expected size AM a ^(x) of the output of TM at b(x) satisfies 

AM a>b (x) ~ , 

1 — x 

so a suitable value of the parameter x to reach target-size n is = 1 — ab/n. 

Lemma 7 (Targeted samplers for the multiset class M). Define 

£ n := 1 - (2C(3)/n) 1 / 3 . 

Then, under the oracle assumption, the expected running time of SampleA^ (£„ ; n) 
is 0(n 4 / 3 ); and, for fixede £ (0, 1), i/ie expected running time o/SAMPLE.A4(£ n ; n j £ ) 
is 0(n 2 / 3 ) as n — > oo, £/ie constant in the big O being independent o/e@. 

In view of stating the expected running times of the targeted samplers for (a x b)- 
boxed multisets, we define the following functions 

0(a) = ^rf. *(a, E ):=^f f (1 + da. 
T{a) T{a) J_ e 

Lemma 8 (Targeted samplers for (a x &)-boxed multisets). Define 

C' b := 1 - ab/n. 

Then, for n > 1, the expected running time of SAMPLEA^ a ,b(Cn' b ! n ) * s equiv- 
alent to 4>(ab)/n as n — > oo. For fixed (a,b,e), the expected running time of 
SAMPLEA4 a .b(£,n' b ; n, e) converges to the constant ao/$(a&, e) as n — » oo. 

3 Precisely, for each fixed e there is no(e) such that the running time is at most en 2 / 3 for 
n > no(e), where the constant c is independent of e. 
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Theorem 9 (Targeted samplers for plane partitions). For n > 1, define the al- 
gorithm S AMPLePartitions [n\ as the procedure that calls S AMPLE.M (£ n ; n) and 
applies Algorithm]^ (Pak's bijection) to the generated diagram. 

Then SAMPLEPARTlTlONS[n] is an exact-size sampler for plane partitions, of 
expected running time 0(n 4 / 3 ). 

For n > 1 and e € (0,1), define SamplePartitions[?t., e] as the algorithm 
that calls Sample.M(£„; n, e) and applies Algorithm [7] to the generated diagram. 
Then SamplePartitions[?t., e] is an approximate- size sampler for plane partitions, 
of expected running time 0(n(lnn) 3 ) as n — > oo (under fixed e), the asymptotic 
constant in the big O not depending on e. 

The proofs of the expected running times announced in Theorem [§] and Theo- 
rem rrU] (given next) are delayed to Section [5] 

In view of stating the expected running times of the targeted samplers for (a x b)- 
boxed plane partitions, we define the following function 

ip(a,b) := \Lt(t + 1) - - 1), where L = max(a, b), £ = min(a,b). 

Theorem 10 (Targeted samplers for (a x 6)-boxed Plane Partitions). For n > 1, 
define SAMPLEPARTlTlONS ai &[n] as the algorithm that calls SAMPLE.M aj &(^' b ; n) 
and applies Algorithm]]^ (Pak's bijection) to the generated diagram. 

Then S AMPLEPARTiTiONS a h [n] is an exact-size sampler for (a x b)-boxed plane 
partitions, of expected running time equivalent to <p(ab)/n as n — > oo. 

Forn > 1 ande G (0, 1), define SAMPLEPARTiTiONS ai f,[n, e] as the algorithm that 
calls SAMPLEA^a, b(Cn' 6 j n i £ ) an ^ applies Algorithm\7\to the generated diagram. 

Then SAMPLEPARTlTlONS aj fc[n, e] is an approximate- size sampler for (a x b)- 
boxed plane partitions, of expected running time equivalent to the constant ip{a, b) + 
ab/Q(ab, e) as n — > oo (under fixed (a, 6, e) ). 

Let us mention that the targeted samplers for (a x &)-boxed plane partitions are 
easily extended to the framework of (a x &)-boxed skew plane partitions. For a fixed 
admissible index-domain D C [0..a — 1] x [Q..b — 1], the appropriate value to reach 
a target size n exactly (or approximately) is ^i £> ' ) := 1 — \D\/n, where \D\ is the 
cardinality of D. 

Another important remark is that the value works well in the asymptotic 
regime, that is, when n » \D\. If not in the asymptotic regime (say one generates 
plane partitions of size 10, 000 constrained to a rectangular box 100 x 100) one 
has to consider the target-size equation more closely. The generating function for 
multisets with support in D is 

n i _ JUi ■ 

Hence the target-size equation — xM' d (x)/Md{x) = n — is 

V (i + j + l)x l +3+ 2 _ 

(hj)eD 

which is to be solved exactly if n is not in the asymptotic regime for the domain D. 
(The solution is asymptotically 1 — \D\/n, but the rate of convergence is slow when 
\D\ is large.) 

5. Analysis of the complexity 



This section is dedicated to proving the expected running times of the random 
samplers, as stated in Theorem |H] and Theorem [TUT Since most of the difficulty is 
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in proving Theorem [5] (unconstrained plane partitions), the proof of Theorem 1 101 is 
only given in the very last subsection (Section 15. 6p and is kept short. 

Recall that the random samplers consist of two steps: generate a diagram under 
the Boltzmann model until the size is in the desired target-domain, and then ap- 
ply Algorithm [T] (Pak's bijection) to the diagram so as to output a random plane 
partition. Accordingly, the complexity of generation is obtained by adding up the 
cost of generating a diagram and the cost of Pak's bijection. 

The expected costs of generating diagrams under Boltzmann model are naturally 
expressed as certain infinite sums, which are best handled by the Mellin transform, 
recalled next. On the other hand, Pak's bijection has complexity cubic in a certain 
parameter called the maximum hook-length, which is the maximal value of the 
hook-length (abscissa+ordinate+1) over all nonzero entries. Therefore, we need to 
find the asymptotic order of the maximum hook-length under the uniform distri- 
bution at size n. 

5.1. The Mellin transform. The Mellin transform is a powerful technique to 
derive asymptotic estimates of expressions involving specific infinite sums, (see [7\ 
for a detailed survey), which occur recurrently in the analysis of our samplers. 
Given a continuous function f(t) defined on R + , the Mellin transform of f(t) is the 
function 



(11) f*(s) := / /(<)i s -Mi. 

Jo 

For instance, the Euler Gamma function T(s) := J °° e _t i s_1 dt is the Mellin trans- 
form of e~*. If f(t) = 0(t- a ) as t -> 0+ and f(t) = 0(^ h ) as t -> +oo, then 
/*(s) is an analytic function defined on the fundamental domain a < Re(s) < b. 
In addition, f*(s) is in most cases continuable to a meromorphic function in the 
whole complex plane (for instance, T(s) is continuable to a meromorphic function 
having its poles at negative integers). In a similar way as the Fourier transform, 
the Mellin transform is almost involutive, the function f(t) being recovered from 
/* (s) using the inversion formula 

/>c+ioo 

(12) /(*)=/ t(s)t- s ds for any cG (a, b). 

J c— zoo 

From the inversion formula and the residue theorem, the asymptotic expansion of 
f(t) as t — > 0~ can be derived from the poles of f*{s) on the left of the fundamen- 
tal domain, the rightmost such pole giving the dominant term of the asymptotic 
expansion. If f*(s) is decreasing very fast as Im(s) — * oo, (which occurs in all the 
series to be analysed next, based on the fact that T(s) is decaying fast and £(s) is 
of moderate growth as Im(s) — > oo), then there holds the following transfer rule [7]: 
a pole of f*(s) of order k+1 (k > 0), 

J ^ s Z a Aa ( s - a) fe +! 

yields a term 

\ a t- a (]nt) k 

in the singular expansion of f(t) around 0. In particular, a simple pole A a /(s — a) 
yields a term \ a /t a . 

Another fundamental property of the Mellin transform is to factorize sums of a 
certain form, 

(13) = 5>*/m s*w = (Ewr)/*w- 

k>l k>l 
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5.2. Complexity of the Boltzmann samplers for multisets. In this section, 
we analyse the complexity of the free Boltzmann sampler TM (x) for the multiset 
class A4, not studying yet the rejection cost when targeting at a certain size-domain. 
In general, given a combinatorial class C for which an explicit Boltzmann sampler 
TC(x) is designed, we write AC(x) for the expected running time of a call to TC(x). 
More generally, we use thereafter the letter A as a prefix to denote the expected 
running time of a random generator. As we are interested in drawing large plane 
partitions, which requires to let x tend to 1, we analyse the asymptotic order of 
AM(x) as x — > 1 — . Recall that A := Z * Seq(Z) 2 , with generating function 
A(x) = x/(l — x) 2 . By definition, the first step of the Boltzmann sampler TM(x) 
is to draw an integer K under the probability distribution 

(14) ¥(K < k) = JJexp(- jA(x j )). 

j>k 

Under the oracle assumption discussed in Section [3] and described in details 
in [5] , the complexity of drawing K is thus of the order of the value k that is finally 
assigned to K. Hence the expected running time of drawing K is of the same order 
as the expected value of K under the above given distribution. 

Lemma 11. The expectation W, X (K) of K under the distribution (I14[) satisfies 

E X (K) = 0((l-xy 1 ln(l-x)) as x -> 1". 

Proof. Fix x G (0, 1). Let r = r(x) be the smallest integer such that x r < (1 — x)/2, 
i.e., r = |hi((l - cc)/2)/ln(x)J + 1. Note that x l < 1/2 for i > r. Hence, for i > r, 
A(x l ) < Ax\ And, for k > r, 

k 

2x k ~ r . 



El ,, , v-^ x % Ax 

7 ^)<4£-< — < 



i>k i~>k 

Hence, for k > r, 

P(K > k) = 1 - P(K <k)<l- cxp(-2x k - r ) < 2x k - r . 
We obtain thus 

2 



E X (K) = F i K >k) <r + 2^ x k - r <r+- 

k>0 k>r 

which concludes the proof since r — r(x) is 0((1 — a;) -1 ln(l — a;)) as x — > 1 _ . □ 

Once the integer K is drawn, the Boltzmann sampler TM(x) draws Poisson 
laws and geometric laws (a call to r^4(a:) consists of two calls to geometric laws). 
Precisely, for each i > 1, the number of calls to TA(x l ) follows a Poisson law 
Pois(A(x i )/i). Since E(Pois(A)) = A, the expected number of calls to TA(x l ) is 
A{x l )/i. In addition, each call to TA takes constant time, since it consists of two 
calls to geometric laws. Hence 

(15) AM(x) = o(E x (K)) +o(^A{x l )h) = +o(^^(a; l )/z). 

i>l i>l 

Lemma 12. The expected running time of the Boltzmann sampler TM(x) satisfies 

AM(x)= O ((1-x)- 2 ). 

Proof. By Equation (|15[) . it is enough to show that F(x) := Y2%>i A(x l )/i is 0((1 — 
x)~ 2 ) as x — » 1~. This is a first instance where the Mellin transform can be 
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successfully applied (more elementary approaches would work in this simple case). 
Define L(t) := F(e~*). Then 

r>l V ' r>l V ; 

The factorization property of the Mellin transform, Equation (|13[) , yields 

r ( s ) = (E; rS )/ , ( s )^( s + 1 )/ , ( s )' 

r>l 

where ((s) :— J2 r >i r S ^ s the Riemann zeta function. Since f(t) = X)n>i ne_ "*' 

the factorization property yields /*(s) = (j2 n>1 nn~ s ^J T(s) = £(s — l)T(s). Thus, 

L*(s) = C(s+l)C(s-l)r(s). It is easily checked that L(t) = 0{l/t 2 ) as t -> 0+ and 
Lit) — 0(e - *) as t — > oo, so that the fundamental domain of L*(s) is Re(s) > 2. 
Hence, to determine the asymptotic behavior of Lit) as t — > + , we have to find 
the rightmost poles of L*(s) such that Re(s) < 2. The function ((s) has a unique 
pole at s — 1 with coefficient 1, and the function T(s) has its poles at non-positive 
integers. Hence L*(s) has a simple pole at s = 2, with coefficient C(3), and no other 
pole for Re(s) > 1, so that the transfer rule of the Mellin transform yields 

Ht) = ^. + o(-) ast^0+. 



t 2 \t, 
The change of variable t = — ln(ar) yields 

Fix) = , ^ 3 \„ + 0( — i— ^ as .t ^ P. 

As a consequence, F(a;) = 0((1 — x)~ 2 ) as x — > l - . □ 

5.3. Analysis of the size of a multiset in M. under the Boltzmann model. 

Given < x < 1, denote by N x the random variable giving the size of the output of 
TM(x) (which is also the size of a plane partition under the Boltzmann model at x); 
Figure H] shows plots of N x for several values of x. As the Boltzmann probability of 
an object of size n is x n /M(x) , the expectation and variance of N x satisfy (see [5] 
for details): 



E(N X ) = Y nM n ^j- = x^-^l, V(N X ) = V n 2 M n ^ T --MN x f .. 

> M ( x ) M(x) ' V ; ^ M(x) V y dx 

Lemma 13. The expectation and variance of the size of a multiset /ieM drawn 
under Boltzmann model satisfy 



(i-z) 3 «-ri- v(i-^) 2 y (i-^) 4 \0--x) 

Proof. We use once again the Mellin transform to derive the asymptotic estimates. 
Observe that M'(x)/M(x) is the logarithmic derivative of M(x), hence the expres- 
sion JTJ of M (x) = P(x) yields 



w = *E^=E'< 



r>l r>l 

Define L(t) :=E(jV e -t). Then 



^)-E- 2 t^ = E- 2 /(^ 
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where f(t) := e-*/(l - e~*) = E n >i e ~"*- Hence 

L*(s) = ^(r 2 r- s )/*( S ) = C(*-2)/* (a) - C(*-2) £ n ^ T ( s ) = C(*-2)C(*)r(*). 
r>l n>l 

The function L*(s) is defined on the fundamental domain Re(s) > 3. The rightmost 
pole such that Re(s) < 3 is at s = 3, where L*(s) ~ 2£(3)/(s — 3). As there are no 
other poles for Re(s) > 2, the transfer rule yields 



L(t) = + 0{r 2 ) as t -> 

Hence the change of variable cc = — ln(£) gives 

The variance is treated similarly, 

dE(i\y 



da; ^ (l-x 7 -) 2 ' 

Hence the function L(t) := Y(N e -t) satisfies L(t) = J2 r >i r3 9( r t)> wnere d(t) = 
e-*/(l - e-*) 2 = E„>i« e "" t - Thus > L *( s ) = C( s ~ 3 )C(s- l)r(s). The location 
of the poles of L* (s) and the transfer rule yields 

L( t) = & + (t- 2 ) ast^0+, 

giving 



□ 



5.4. Complexity of the targeted samplers for multisets. Recall that the 
targeted samplers for the multiset class M. repeat calling the Boltzmann sampler 
TM (x) with a suitable value of x until the size is in the target domain f2; SI = {n} 
for exact-size sampling and Q = [n(l — e), n(l + s)] for approximate-size sampling. 

Lemma 14. For n > 1. /e< 6e i/ie solution o/2£(3)/(l — a;) 3 = n, i.e., 

€„ = 1 - (2C(3)/n) 1 /3. 

Define tt u as the probability that the output o/I\M(£ n ) has size n. For any e € (0, 1), 
define i^n,e as the probability that the size of the output afTM(£ n ) is in the range 
[n(l — e),7i(l +e)]. Then, tt„ ~ Cn~ 2 / 3 as n — > oo, with C ~ 0.1082; and, for fixed 
e £ (0, 1), 7r„ i£ — ► 1 as n — > oo. 

Proof. As £ n is solution of 2£(3)/(l — a:) 3 = n, Lemma [TBI ensures that E(A^ n ) = 
n + 0{n 2 l z ) as in oo, i.e., there exists C > such that |E(iVf n ) — n\ < Cn 2 / 3 . 
Hence Chebyshev's inequality gives, for any e G (0, 1), 

l-7r„, e = F(\N U - n\ > en) 



< 



(| W{ .- E(A r { .)|> (e „- C „ 2 /»))<^<gJ 



(en — Cn 2 / 3 ) 2 

Given the fact that V(%J = 0((1 - ^„) -4 ) = 0(n 4/3 ), we have 1 - tt„ >£ -► as 
n — » oo. 

Next we prove the estimate of tt„. Note that 

rr„ = M n ■ (£n) n /M(Z n ) = P n ■ (e„) n /P(C„). 

Hence it is enough to obtain the asymptotics of P n and (£ n ) n as n — > oo, and of 
P{x) as a; — > l - . These have first been found by Wright [H] and later by Meinardus 
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in a more general framework [T7] relying on the saddle-point method (a detailed 
and accessible presentation of the saddle-point method is given in [5], Ch.VIII], 
partitions are studied in the 6th section of the chapter). We briefly review the 
main ingredients. To find the asymptotics of P(x) as x — > 1 — , one applies the 
Mellin transform techniques to the series L(t) = ln(P(e~')), and finds 



(Hi) P(x)^C 1 (l-x) i / 12 exp(^(3)-^ r ^) ms, -I 



with d an explicit constant, d w 0.9368. Define c := (2C(3)) 1 / 3 , so C(3) = c 3 /2 
and £ n = 1 — cnT 1 ! 3 . From we obtain 

Pfe) ~ Cin-^exp (±cn 2 / 3 - ±cV/ 3 ), 

with C[ = Cic 1 / 12 Ri 0.9599. The asymptotics of (£ n ) n is easy to obtain. Since 
log(l — cn -1 / 3 ) = — cn -1 / 3 — ic 2 ?! -2 / 3 — ^-n -1 + o(n _1 ), we obtain 

(£„)" = exp(nlog(l - cn" 1 / 3 )) - C 2 exp ( - cn 2 / 3 - ±cV/ 3 ), 

with C2 = exp(— c 3 /3) « 0.4487. Finally the asymptotics of P n has been obtained 
by Wright [29] using the saddle-point method and the estimate (fT6| . The idea is 
to use Cauchy's formula 



Pn = ^~ I P{z) Z - n - 1 &Z 1 

/c(o,e„) 



2i?r 



with the circle of radius £„ centered at as the integration contour. Using the 
estimate (|16p (more precisely one needs the fact that this estimate holds in an 
open cone centered at 1 and containing the line z < 1), one shows that the main 
contribution of the integral is on a small arc of C(0,£„) around the origin, and 
obtains 

P n ~ C 3 n- 25 / 36 exp(|cn 2/3 ), 
with C3 ?» 0.2315. Thus, from the estimates of P(£n), (£n) n , and P„, we find 

7T„ - C71- 2 / 3 , 

with C w 0.1082. □ 

It is easily checked that the expected running time of a rejection sampler is the 
expected running time of the sampler times the expected number of calls (which is 
the inverse of the probability of success), therefore 

AA/ (£„) A ^ c __ A „ r , . n A AM(£ n ) 



A(sAMPLEA4[C„;n]) = — ^22, a(SampleA<[6 



Since AAf (x) = 0((1 - x)~ 2 ) and 1 - £„ = O^- 1 / 3 ), we have AM(£„) = 0(n 2 / 3 ). 
Moreover — 0(n 2 / 3 ) and l/ir n ,e —* 1 (for fixed e S (0, 1)) bv LemmafT4l Hence 
A(Sample.M(£„; n)) = 0(n 4 / 3 ) and A(SampleA1(^„; n, e)) = <3(n 2 / 3 ), which con- 
cludes the proof of Lemma [7j 

5.5. Complexity of Pak's bijection. The aim of this section is to provide an O 
bound on the expected running time of Algorithm [T] for a multiset /1 e M of 
size n taken uniformly at random. As we will see, the complexity of Algorithm [T] 
applied to a multiset is expressed in terms of the width and length of the bounding 
rectangle of \i. These parameters have been recently studied by Mutafchiev in (TB] : 
using the saddle-point method he shows that the width (similarly the length) after 
suitable normalization converges weakly to an explicit distribution. Since we are 
only interested in a big O, we only need upper bounds, which are much simpler to 
show. Therefore we prefer to provide here our own simple self-contained analysis. 
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For a multiset (j, € M := MSet(Z x Seq(Z) 2 ) represented by its diagram, let w 
and h be the width and height of the bounding rectangle of /i. Pak's bijection 
scans the double range [0 < i < w — 1,0 < j < h — 1]; when a square is 
treated, the squares that are updated are those on the up-right diagonal + c, j + 
c) such that i + c < w, j + c < h}; each update of an entry consists of a fixed number 
of operations involving {+, — , max, min}. The sum of the lengths of the up-right 
diagonals over the squares of the bounding rectangle is X)S!i i(w — i + h—i+1), 
which is equal to 

(17) ijj(w, h) := \LI{1 + 1) - i(^ 3 - 1), where L = max(w, h), i = min(w, h). 

This quantity is clearly 0(L 3 ), so that the complexity of Algorithm [T] is cubic in L. 

We introduce a parameter that will crucially simplify the analysis. Given /i 6 
M. represented as a diagram, the hook-length of an entry of the diagram is 
defined as h(i,j) := i + j + 1, i.e., h(i,j) is the size of seen as an element 
of A = Z x Seq(.Z) 2 . The maximum hook-length of /i, denoted by is the 

maximum value of the hook-length over all non-zero entries of the diagram of /x. 

Lemma 15. Given an element /i in M := MSet(Z x Seq(Z) 2 ), the complexity 
of Pak's bijection applied to /i is 0([k(p)} 3 ), where fc(/i) is the maximal hook-length 
of ii. 

Proof. The maximal hook-length fc(/i) is at least equal to the width w and to the 
height h of the bounding rectangle of [i. Hence the complexity of Pak's bijection, 
which is 0([max(w, h)} 3 ), is also 0([/c(^)] 3 ). □ 

Hence, to have a big O bound on the expected running time of Pak's bijection, 
we need to bound the expected value of k(fi) 3 for a multiset [i £ M. of size n 
taken uniformly at random. First, for a series C(x) — ^2 n c n x n with non-negative 
coefficients and with radius of convergence p > 0, we recall the trivial bound 

(18) c n <C{x)x~ n for any x€ (0,p). 

Lemma 16 (expected running time of Pak's bijection at a fixed size). For n > 1, 

let p e M be a multiset of size n taken uniformly at random. Then the expected 
running time of Algorithm^ applied to /i is 0(n(lnn) 3 ). 

Proof. Denote by H n the expectation of k(/j,) 3 for a multiset fi 6 A4 of size n taken 
uniformly at random. By Lemma 1151 the expected running time of Algorithm [1] 
under the uniform distribution (on M.) at size n is 0(H n ). Hence to show the 
lemma we just have to show that H n = 0(n(lnn) 3 ). For k > 1, denote by 
the family of multisets in Ai with maximal hook-length equal to fc, and denote by 
AfW(x) the series of M {k) . Define 

K{x) :=J]fc 3 M (fc) (x). 

k>l 

Note that H n = [x n ]K(x)/[x n ]M{x) = [x n ]K (x)/P n , with P n the number of plane 
partitions of size n. For k > 1, define a k-pointed multiset as a multiset /i 6 M 
where a non-zero entry at hook-length A; is marked in the diagram of fi. Note that 
the series of fc-pointed multisets is kx k M(x), where the factor k counts the possible 
places to mark an entry and where the factor x k takes account of the fact that the 
marked entry is non-zero. Since fc-pointed multisets form a superfamily of MS k \ 
we have 

M {k \x) < kx k M{x) = kx k P(x), for any < x < 1. 
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Let B be a constant whose value is to be fixed later, and define u n := [Bii 1 / 3 log(n)J . 
We have 

i n 3 71 

H n = ^fc 3 ^MW(x) < K) 3 + \x n ]M {h \x). 

k—1 k—u n 

As M^ k \x) < kx k P(x), the bound (fT5j) ensures that, for u n < k < n, 

[^]M«(x) < M (fc) (Cn)-(^)" n < A;(en) fc -P(en)-(en)" n < n(k)*» • P(£ n ) " (£n)~* 

Hence 

gn<K) 3 +n 5 (^)"" p P( f;\ n . 
Let c := (2C(3)) 1 / 3 . We have 

(&,)"" = exp (u„log(l - en" 1 / 3 )) ~ exp(-cu„n- 1 / 3 ) - n~ cB . 
Moreover, according to Lemma fTH 



= 0{n z ' 3 ). 



Bn ' (<^n ) n TTn 

Hence 

ff„ = 0(n(mn) 3 ) + 0(n 5 n~ ci3 n 2/3 ). 

Taking the constant B sufficiently large so that cB > 4 + 2/3 (e.g., B — 5), we 
obtain 

H n = 0(n(lnn) 3 ). 

□ 

Proposition 17. For any e e (0, 1), i/ie expected running time o/SAMPLEPARTiTiONS[n, e] 
satisfies 

A^SamplePartitions[tj, e]j = 0(n(lnn) 3 ) as n — > oo, 

t/ie asymptotic constant in the big O being independent of e. 
The expected running time of SamplePartitions^] satisfies 

A amplePartitions [ti] J = 0(n 4 / 3 ) as n — > oo. 

Proof. Start with the proof for the exact-size sampler. By definition, we have 

A ( S amplePartitions [n] j = A( SampleA4[^„; u] ) + E„(PakBijection), 



where E„(PakBijection) is the expected running time of Algorithm [T] (Pak's bi- 
jection) for a multiset of size n taken uniformly at random. Lemma [7] ensures that 
A(SAMPLE.M[£ n ;ra]) is 0(n 4 / 3 ); by Lemma [HA E„(PakBijection is 0(n(\nn) 3 ). 
Hence A(SAMPLEPARTlTlONS[n]) is 0(n 4 / 3 ). 

Consider now the approximate-size sampler. By definition, we have 

A^SAMPLEPARTiTiONs[n,e]^ = A ^SampleX [£„ ; n, e]J + E„ !£ (PakBijection), 

where E nj£ (PAKBiJECTiON) is the expected running time of Pak's bijection for a 
multiset drawn from SAMPLE.M(£ n ; n, e). By Lemma[71 A(SAMPLEj\4(£ n ; n, e)) is 
0(n 2//3 ). Since the size of an object output by SAMPLE.M (£ n ; n, e) is at most 2n 
(because e 6 (0,1)), Lemma fTrJl ensures that E njE (PakBijection) is 0(n(hm) 3 ). 
Hence A(SamplePartitions[ti, e]) is 0(n(lnn) 3 ). □ 
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5.6. Complexity of the samplers for (a x &)-boxed plane partitions. By 

definition (see Algorithm [5|) , the Boltzmann sampler TM a ^(x) for (a x 6)-boxed 
multisets just consists of ab calls to geometric laws. Giving unit cost to a call to a 
geometric law, one has 

AM a , b (x) = ab. 

Next, recall that 

n^i+j + l c ^ 
■ , ■ i i ~ -, r~r, with c = I . 
1 - x l +3+ l x ^i- 1 - x) ab 11 i + j + 1 

0<i<a v ' 0<i<a J 

0<j<b 0<j<b 

For a > 0, a class C is called a-singular at x = 1 if C(x) ~ c/(l — a;)" for some 
constant c > (the ~ holding in a complex neighbourhood of 1) and if 1 is the only 
singularity of C(x) in a disk of the form {zeC s.t. |z| < 1 + S} for some (5 > 0. 
Note that M a ,b is a-singular for a — ab. In [5] it is shown that if C is a-singular 
at x = 1 , then for n > 1 the probability 7r„ of being of size n under the Boltzmann 
model at x n := 1 — a/n satisfies 

0(a) , ,/ \ (a/e) a 
(19) 7r„ ~ , where 0(a) - 



a— 1 „-qs. 



And, for fixed e £ (0, 1), the probability 7r„. e of being in the size-domain [n(l 
e),n(l + e)] under the Boltzmann model at x n satisfies 

(a/e) a f E 

(20) 7r n , E -» $( a ,e), where $(a, e) = i-f^- / (1 + s) 

Since the expected running time of a rejection sampler is the expected running time 
of the sampler divided by the probability of success at each attempt, the expected 
running times of the targeted samplers for M a ,b satisfy asymptotically 

A(SAMPLE.M a , 6 [C & ;™]) ~ TTIT^ A(SAMPLEX a , h [C'>,£]) -► ^ v 
\ J n-*oo (p[aO) V / n^oo q>(ab, £) 

The second step of the targeted samplers for (a x 6)-boxed plane partitions is 
Algorithm [T] (Pak's bijection). When x — > 1~, all entries of the rectangle i? aj f, := 
[0..a — 1] x [0..b — 1] in the diagram of fj, <— TM a ^{x) are non-zero with high 
probability, hence the bounding rectangle of /i is R a ,b with high probability. As a 
consequence, the complexity of Algorithm [T] applied to /i is with high probability 
the quantity ip(a,b) defined in (jTTJ) . Hence 

A TSAMPLEPARTITIONSa.b [n] \ ~ °^ 11 + lp(a, b) 



4>(ab) ' n-^oo (j){ab) 



A( SamplePartitionSq &[n,e] — * — — — -+ip(a,b), 
\ ' / n— >oo q>(ao, e) 

which concludes the proof of Theorem fTUl 
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